Se requieren cargar las siguientes paqueterías para seguir los ejemplos.
Los modelos que veremos aquí tienen como idea principal encontrar relaciones lineales entre la serie que queremos pronosticar, \(y\), con una o más series distintas, x. En otras palabras, pronosticaremos los valores futuros de una serie, a partir de los cambios en otra serie que la afecte.
Es muy común querer predecir de esta forma. Por ejemplo, una tienda de helados podría encontrar una relación entre sus ventas (\(y\)) y la temperatura (\(x_1\)). O las ventas de Nike, a partir de cuánto gastan en publicidad y mercadotecnia.
En la literatura podemos encontrar muchos nombres para las variables \(y\) ^ \(x\). P. ej.
| \(y\) (var. de pronóstico) | \(x\) (vars. predictoras) |
|---|---|
| Var. dependiente | Vars. independientes |
| Explicada | Explicativas |
| Regresada | Regresoras |
| Respuesta | Estímulo |
| Resultado | Covariante |
| Controlada | De control |
El caso más sencillo sería un modelo de regresión lineal simple, de la forma:
\[ y_t = \beta_0 + \beta_1 x_t + \varepsilon_t \]
donde
\(\beta_0\) es conocido como el intercepto y representa el valor predicho cuando \(x = 0\).
\(\beta_1\) es la pendiente de la recta. Nos indica el cambio promedio en \(y\), ante un cambio en una unidad de \(x\).
El término de error, \(\varepsilon_t\) se asume aleatorio y decimos que captura los cambios debido a todas las otras variables que pudieran llegar a afectar a \(y_t\), que no están explícitamente especificadas en el modelo.
La recta resultante está dada entonces por \(\beta_0 + \beta_1 x_t\), y la diferencia que existe en los puntos reales y ésta es \(\varepsilon_t\).
Como primer ejemplo, veamos las tasas de crecimiento del gasto de consumo, \(y\), y su relación con el ingreso personal disponible, \(x\).
La gráfica de tiempo de ambas series:
us_change %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Consumption, colour = "Consumo")) +
geom_line(aes(y = Income, colour = "Ingreso")) +
ylab("cambio %") + xlab("Año") +
guides(colour=guide_legend(title="Series")) +
theme(legend.position = "top")Un diagrama de dispersión entre ambas series, para ver una posible correlación.
us_change %>%
ggplot(aes(x=Income, y=Consumption)) +
ylab("Consumo (cambio % trimestral)") +
xlab("Ingreso (cambio % trimestral)") +
geom_point() +
geom_smooth(method="lm", se=FALSE)La línea azul en el gráfico sigue la ecuación que describe la regresión lineal que tiene por variable dependiente al consumo y como independiente al ingreso.
El término de regresión fue acuñado por primera vez por Francis Galton en 1886. Él estaba estudiando la relación que existe entre la estatura de los hijos y la estatura de los padres.
Lo que encontró fue lo siguiente, en resumen:
Los padres más altos, tendían a tener hijos más altos, mientras que los padres bajos tendían a tener hijos bajos.
En promedio, los hijos de padres altos no logran ser más altos que ellos. Similarmente, los hijos de padres bajos, en promedio son más altos que sus papás.
Así, Galton decía que había una tendencia a regresar a la estatura promedio.
De no cumplirse la regresión de Galton, sería común tener gente de la estatura de un Hobbit, y también de la estatura de un gigante.
Entonces, el análisis de regresión en tiempos modernos trata sobre la relación de la dependencia entre una variable \(y\), respecto de una o más variables exógenas (regresoras \(x\)) para predecir el valor promedio de la variable dependiente.
“Una relación estadística, por más fuerte y sugerente que sea, nunca podrá establecer una conexión causal: nuestras ideas de causalidad deben provenir de estadísticas externas y, en último término, de una u otra teoría” (Kendall & Stuart, 1961)
Regresión \(\neq\) Causalidad
Correlación \(\neq\) Causalidad
Se puede hablar de linealidad en dos sentidos:
Para un modelo de regresión lineal, solo nos interesa que sea lineal en los parámetros.
Todas estas funciones son lineales en los parámetros y pueden ser estimadas mediante un modelo de regresión lineal
Así, un modelo de regresión lineal puede generar una recta, o una variedad de curvas, dependiendo la forma funcional que se elija.
\[ H_0: \beta_0 = 0 \\ H_1: \beta_0 \neq 0 \] \[ H_0: \beta_1 = 0 \\ H_1: \beta_1 \neq 0 \]
\[ \hat{y}_t = \hat{\beta}_0 + \hat{\beta}_1x_t + \hat{\varepsilon}_t \\ \hat{y}_t = 0.54454 + 0.27183 x_t + \hat{\varepsilon}_t \]
\[ y_{consumo} = \beta_0 + \beta_1 x_{income} + \beta_2 x_{production} + \beta_3 x_{savings} + \beta_4 x_{unemployment} \]
Una correlación, por más fuerte que sea entre dos variables, no puede implicar por sí misma causalidad.
us_change %>%
pivot_longer(cols = -Quarter) %>%
ggplot(aes(x = Quarter, y = value, color = name)) +
geom_line() +
facet_wrap(~ name, scales = "free_y") +
theme(legend.position = "none")Graficando nuestra variable de pronóstico (Consumo) vs. cada una de las variables predictoras:
us_change %>%
pivot_longer(cols = -c(Quarter, Consumption)) %>%
ggplot(aes(x = Quarter, y = value, color = name)) +
geom_line() +
geom_line(aes(y = Consumption), color = "black") +
facet_wrap(~ name, scales = "free_y") +
theme(legend.position = "none")Realizamos un primer modelo, donde utilizaremos de variable predictora al ingreso disponible, para pronosticar el consumo.
\[ y_{t,Consumo} = \beta_0 + \beta_{1}x_{t,Ingreso} + \varepsilon_t \]
Series: Consumption
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-2.58236 -0.27777 0.01862 0.32330 1.42229
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.54454 0.05403 10.079 < 2e-16 ***
Income 0.27183 0.04673 5.817 2.4e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5905 on 196 degrees of freedom
Multiple R-squared: 0.1472, Adjusted R-squared: 0.1429
F-statistic: 33.84 on 1 and 196 DF, p-value: 2.4022e-08
Prueba de significancia individual para cada uno de los parámetros (\(\beta_i\))
\[ H_0: \beta_i = 0 \]
Prueba de significancia conjunta (para revisar si nuestro modelo sirve o no)
\[ H_0: \beta_1 = \beta_2 = \beta_3 = \ldots = 0 \]
augment(fit1) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Consumption, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted"))+
xlab("Año") + ylab(NULL) +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
guides(color = guide_legend(title = NULL))El modelo no parece capturar adecuadamente la variación de los datos reales.
augment(fit1) %>%
ggplot(aes(x = Consumption, y = .fitted)) +
geom_point() +
ylab("Fitted (valores ajustados)") +
xlab("Datos (reales históricos)") +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
geom_abline(intercept = 0, slope = 1)Es evidente que este modelo se puede mejorar. Probemos incluyendo las otras predictoras.
Podríamos proponer ahora un modelo en donde el consumo sea una función del ingreso, la producción, los ahorros y el desempleo:
\[ y_{t,Consumo} = \beta_0 + \beta_{1}x_{t,Ingreso} + \beta_2x_{t,Producción} + \beta_3x_{t,Ahorros} + \beta_4x_{t,Desempleo} + \varepsilon_t \]
fit2 <- us_change %>%
model(
reg_lin_multiple = TSLM(Consumption ~ Income + Production + Savings + Unemployment)
)
report(fit2)Series: Consumption
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-0.90555 -0.15821 -0.03608 0.13618 1.15471
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.253105 0.034470 7.343 5.71e-12 ***
Income 0.740583 0.040115 18.461 < 2e-16 ***
Production 0.047173 0.023142 2.038 0.0429 *
Savings -0.052890 0.002924 -18.088 < 2e-16 ***
Unemployment -0.174685 0.095511 -1.829 0.0689 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3102 on 193 degrees of freedom
Multiple R-squared: 0.7683, Adjusted R-squared: 0.7635
F-statistic: 160 on 4 and 193 DF, p-value: < 2.22e-16
augment(fit2) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Consumption, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted"))+
xlab("Año") + ylab(NULL) +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
guides(color = guide_legend(title = NULL))Este modelo parece capturar más variación de los datos históricos.
augment(fit2) %>%
ggplot(aes(x = Consumption, y = .fitted)) +
geom_point() +
ylab("Fitted (valores ajustados)") +
xlab("Datos (reales históricos)") +
ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
geom_abline(intercept = 0, slope = 1)df <- left_join(us_change, residuals(fit2), by = "Quarter")
df %>%
select(-c(Consumption, .model)) %>%
pivot_longer(cols = c(Income:Unemployment)) %>%
ggplot(aes( x = value, y = .resid, color = name)) +
geom_point() + ylab("Residuales") + xlab("Predictoras") +
facet_wrap(~ name, scales = "free_x") +
theme(legend.position = "none")augment(fit2) %>%
ggplot(aes(x = .fitted, y = .resid)) +
geom_point() +
labs(x = "Ajustados", y = "Residuales")fit3 <- us_change %>%
model(r1 = TSLM(Consumption ~ Income),
r2 = TSLM(Consumption ~ Income + Production),
r3 = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
r4 = TSLM(Consumption ~ Income + Production + Savings),
r5 = TSLM(Consumption ~ Income + Savings + Unemployment),
r6 = TSLM(Consumption ~ Income + Production + Unemployment),
r7 = TSLM(Consumption ~ Income + Savings)
)
fit3 %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)Series: Consumption
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-0.90555 -0.15821 -0.03608 0.13618 1.15471
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.253105 0.034470 7.343 5.71e-12 ***
Income 0.740583 0.040115 18.461 < 2e-16 ***
Production 0.047173 0.023142 2.038 0.0429 *
Savings -0.052890 0.002924 -18.088 < 2e-16 ***
Unemployment -0.174685 0.095511 -1.829 0.0689 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3102 on 193 degrees of freedom
Multiple R-squared: 0.7683, Adjusted R-squared: 0.7635
F-statistic: 160 on 4 and 193 DF, p-value: < 2.22e-16
Selección del modelo a través de backwards stepwise regression con base en el \(AICc\) sugiere conservar el modelo que incluye todas las predictoras.
us_change %>%
model(i = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
ii = TSLM(Consumption ~ Income + Production + Savings),
iii = TSLM(Consumption ~ Income + Production + Unemployment),
iv = TSLM(Consumption ~ Income + Savings + Unemployment),
v = TSLM(Consumption ~ Production + Savings + Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)Con base en \(BIC\) parece que el mejor modelo es el que incluye todas las predictoras menos el desempleo.
us_change %>%
model(i = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
ii = TSLM(Consumption ~ Income + Production + Savings),
iii = TSLM(Consumption ~ Income + Production),
iv = TSLM(Consumption ~ Income + Savings),
v = TSLM(Consumption ~ Production + Savings),
vi = TSLM(Consumption ~ Income + Savings + Unemployment),
vii = TSLM(Consumption ~ Production + Savings + Unemployment),
viii = TSLM(Consumption ~ Income + Production + Unemployment)
) %>%
glance() %>%
select(.model, adj_r_squared, AIC, AICc, BIC)En estos pronósticos solo se utiliza información disponible hasta el último dato del histórico. A estos pronósticos se les considera como pronósticos reales. Aquí las predictoras se deben pronosticar antes de poder producir el pronóstico de la variable de interés.
Con estos pronósticos se utiliza información real disponible de las predictoras. Estos pronósticos ya no son reales (en el sentido estricto). La variable a pronosticar (\(y\)) sigue siendo desconocida.
fit_escenarios <- us_change %>%
model(lineal = TSLM(Consumption ~ Income + Savings + Unemployment))
# Necesitamos agregar nuevos datos de las predictoras
optimista_futuro <- new_data(us_change,4) %>%
mutate(Income = 0.5,
Savings = 0.5,
Unemployment = 0)
pesimista_futuro <- new_data(us_change,4) %>%
mutate(Income = -1,
Savings = -0.5,
Unemployment = 0)
fc_optimista <- forecast(fit_escenarios, new_data = optimista_futuro) %>%
mutate(Escenario = "Optimista") #%>% # REVISAR
# as_fable(response = "Consumo", distribution = Consumption, key = c("Escenario",".model"))
fc_pesimista <- forecast(fit_escenarios, new_data = pesimista_futuro) %>%
mutate(Escenario = "Pesimista") # %>% # REVISAR
# as_fable(response = "Consumo", distribution = Consumption, key = c("Escenario",".model"))
# bind_rows(fc_optimista,fc_pesimista)
# us_change %>%
# autoplot(Consumption) +
# autolayer(bind_rows(fc_optimista,fc_pesimista))
p1 <- us_change %>%
autoplot(Consumption) +
autolayer(fc_optimista) +
ggtitle("Escenario optimista")
p2 <- us_change %>%
autoplot(Consumption) +
autolayer(fc_pesimista) +
ggtitle("Escenario pesimista")
p1 / p2Vamos a generar un pronóstico ex-ante del consumo. Para esto, necesitamos primero producir pronósticos de las predictoras:
# mod_predictoras <- function(predictora, horizonte = 4) {
# us_change %>%
# model(predictora = ARIMA(as.formula(predictora)) %>%
# forecast(h = horizonte)
# }
#
# mod_predictoras(predictora = Income)
ingreso <- us_change %>%
model(ETS = ETS(Income),
ARIMA = ARIMA(Income)
) %>%
forecast(h = 4)
ingreso %>%
autoplot(us_change, level = NULL)
produccion <- us_change %>%
model(
ETS = ETS(Production),
ARIMA = ARIMA(Production)
) %>%
forecast(h = 4)
produccion %>%
autoplot(us_change, level = NULL)
ahorro <- us_change %>%
model(
ETS = ETS(Savings),
ARIMA = ARIMA(Savings)
) %>%
forecast(h = 4)
ahorro %>%
autoplot(us_change, level = NULL)
desempleo <- us_change %>%
model(
ETS = ETS(Unemployment),
ARIMA = ARIMA(Unemployment)
) %>%
forecast(h = 4)
desempleo %>%
autoplot(us_change, level = NULL)Teniendo ya los pronósticos de cada predictora, podemos proceder a generar el pronóstico del consumo.
fit <- us_change %>%
model(
`Regresión lineal múltiple` = TSLM(Consumption ~ Income + Production + Savings + Unemployment)
)
datos_futuros <- new_data(us_change,4) %>%
mutate(Income = ingreso %>% filter(.model == "ARIMA") %>% pull(.mean),
Savings = ahorro %>% filter(.model == "ARIMA") %>% pull(.mean),
Unemployment = desempleo %>% filter(.model == "ARIMA") %>% pull(.mean),
Production = produccion %>% filter(.model == "ARIMA") %>% pull(.mean))
datos_futuros
recent_production %>%
autoplot(Beer) +
labs(x = "Año", y = "Megalitros",
title = "Producción de cerveza trimestral en Australia")Existen varias predictoras que pueden ser útiles en el análisis de regresión.
trend()\[ y_t = \beta_0 + \beta_1t + \varepsilon_t \]
Las variables dummy (también llamadas dicotómicas, binarias, …) toman solo dos valores: 1 o 0 (verdadero/falso).
Para agregar variables dummy estacionales, basta con escribir season().
recent_production %>%
select(Quarter,Beer) %>%
mutate(tendencia = seq_along(recent_production$Quarter),
q2 = if_else(quarter(Quarter)==2,1,0),
q3 = if_else(quarter(Quarter)==3,1,0),
q4 = if_else(quarter(Quarter)==4,1,0)
) %>%
model(TSLM(Beer ~ tendencia + q2 + q3 + q4)) %>%
report()Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-42.9029 -7.5995 -0.4594 7.9908 21.7895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 441.80044 3.73353 118.333 < 2e-16 ***
tendencia -0.34027 0.06657 -5.111 2.73e-06 ***
q2 -34.65973 3.96832 -8.734 9.10e-13 ***
q3 -17.82164 4.02249 -4.430 3.45e-05 ***
q4 72.79641 4.02305 18.095 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12.23 on 69 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9199
F-statistic: 210.7 on 4 and 69 DF, p-value: < 2.22e-16
Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-42.9029 -7.5995 -0.4594 7.9908 21.7895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 441.80044 3.73353 118.333 < 2e-16 ***
trend() -0.34027 0.06657 -5.111 2.73e-06 ***
season()year2 -34.65973 3.96832 -8.734 9.10e-13 ***
season()year3 -17.82164 4.02249 -4.430 3.45e-05 ***
season()year4 72.79641 4.02305 18.095 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12.23 on 69 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9199
F-statistic: 210.7 on 4 and 69 DF, p-value: < 2.22e-16
p <- augment(fit_beer) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Beer, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted")) +
labs(x = "Año", y = "Megalitros")
ggplotly(p)Al ver la gráfica, vemos que existe un trimestre muy por debajo del resto. Podemos agregar una variable de intervención para ese periodo.
Capturan el efecto de un solo periodo.
cerveza <- recent_production %>%
select(Quarter, Beer) %>%
mutate(q2_94 = if_else(Quarter == yearquarter("1994 Q2"),1,0),
q4_92 = if_else(Quarter == yearquarter("1992 Q4"),1,0))
cervezaAjustamos un modelo, corrigiendo por ese periodo outlier (1994 Q2).
Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-42.6171 -6.9644 -0.8231 8.1306 21.7895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 442.55585 3.67748 120.342 < 2e-16 ***
trend() -0.36068 0.06603 -5.462 7.20e-07 ***
season()year2 -33.34433 3.94421 -8.454 3.28e-12 ***
season()year3 -17.82164 3.94060 -4.523 2.51e-05 ***
season()year4 72.81682 3.94115 18.476 < 2e-16 ***
q2_94 -24.60467 12.46254 -1.974 0.0524 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 11.98 on 68 degrees of freedom
Multiple R-squared: 0.9284, Adjusted R-squared: 0.9232
F-statistic: 176.4 on 5 and 68 DF, p-value: < 2.22e-16
Capturan el efecto a partir de cierto periodo.
Series: Beer
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-43.0235 -7.5945 -0.4539 7.7754 21.4829
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 441.9154 3.8644 114.354 < 2e-16 ***
trend() -0.3548 0.1308 -2.712 0.00846 **
season()year2 -34.6452 3.9985 -8.665 1.36e-12 ***
season()year3 -17.8046 4.0536 -4.392 4.02e-05 ***
season()year4 72.8279 4.0594 17.941 < 2e-16 ***
d2000 0.7282 5.6402 0.129 0.89766
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12.32 on 68 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9188
F-statistic: 166.1 on 5 and 68 DF, p-value: < 2.22e-16
p <- augment(fit_beer) %>%
ggplot(aes(x = Quarter)) +
geom_line(aes(y = Beer, color = "Datos")) +
geom_line(aes(y = .fitted, color = "Fitted")) +
labs(x = "Año", y = "Megalitros")
ggplotly(p)Cuando se crean varaibles dummy, la cantidad de dummies a generar es una menos que el total de categorías. Si son dos categorías, solo necesitamos una dummy. Si son tres, necesitamos 2, etc.
\[ \log{y_t} = \beta_0 + \beta_1 \log{x_t} \]
La interpretación de los coeficientes (\(\beta_s\)) es como elasticidades (cambios porcentuales).
\[ y_t = \beta_0 + \beta_1 \log{x_t} \]
\[ \log{y_t} = \beta_0 + \beta_1 x_t \]
Aquí la regresión se calcula por partes en el tiempo, para capturar distintas tendencias.
boston_men <- boston_marathon %>%
filter(Event == "Men's open division") %>%
mutate(Minutes = as.numeric(Time)/60)
boston_men %>%
autoplot(Minutes) +
ggtitle("Tiempos ganadores del maratón de Boston, categoría abierta de hombres")Vamos a modelar los tiempos con una regresión por partes.
fit_boston <- boston_men %>%
model(
lineal = TSLM(Minutes ~ trend()),
exponencial = TSLM(log(Minutes) ~ trend()),
`Reg. por partes` = TSLM(Minutes ~ trend(knots = c(1940,1980)))
)
fc_boston <- fit_boston %>% forecast(h = 10)
boston_men %>%
autoplot(Minutes) +
geom_line(aes(y = .fitted, color = .model), data = fitted(fit_boston)) +
autolayer(fc_boston, alpha = 0.5, level = 95) +
ggtitle("Maratón de Boston, cat. abierta de hombres")